% albedos = getAlbedos(lights, normals, scImgs, scMask)
%
% input:
%   * lights - A 3-by-numImages matrix representing the lighting direction
%   of each image in the image series. lights(:,i) is the lighting
%   direction vector of the ith image in the series. A lighting direction
%   vector is of the form (l_x, l_y, l_z), where the x axis is along the
%   column index increasing direction (towards right) and the y axis is
%   along the row index increasing direction (downwards). The z axis points
%   out of the image.
%   * normals - A numRows-by-numColumns-by-3 matrix representing the
%   normals of each pixel in a numRows-by-numColumns image. normals(i,j,:)
%   is the normal vector of the scene point corresponding to pixel(i,j) in
%   the image. The normal vector is in the form of (n_x, n_y, n_z), where
%   the x axis is along the column index increasing direction (towards
%   right) and the y axis is along the row index increasing direction
%   (downwards). The z axis points out of the image.
%   * scImgs - A numRows-by-numColumns-by-numImages matrix representing a
%   series of images taken for the same scene under different lighting
%   conditions. scImgs(:,:,i) is the ith image in the image series, 
%   presented in a single color channel.
%   * scMask - A numRows-by-numColumns matrix masking the scene related
%   pixels. However, this matrix doesn't consist of only 0 and 1 indicating
%   exactly whether a pixel is or is not a scene point. The values kind of
%   range from 0 to 255. In this method, if scMask(i,j) is larger than
%   127.5, then pixel(i,j) in the image is considered to be a scene point;
%   otherwise, the background pixel.
%
% output:
%   * albedo - A numRows-by-numColumns matrix representing the albedo for
%   each scene point in the image (under single color channel).
%
% description:
% This method takes the lighting directions, scene point normals, a
% single-color image series under different lighting conditions, and a
% single-color image mask of scene points to calculate the albedos for each
% scene point pixel.

function albedos = getAlbedos(lights, normals, scImgs, scMask)
[nRs, nCs, nImgs]  = size(scImgs);
albedos = zeros(nRs, nCs);
threshold = max(max(scMask))/2;
for i = 1:nRs
    for j = 1:nCs
        if scMask(i,j) >= threshold
            n = reshape(normals(i,j,:),3,1);
            intensities = double(reshape(scImgs(i,j,:),1,nImgs));
            J = -lights'*n;
            num = intensities*J;
            den = J'*J;
            k = num/den;
            albedos(i,j) = k;
        end
    end
end
end

